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ABSTRACT 


A  sixteen  node  shell  element  is  developed  using  a  matrix  stabilization 
scheme  based  on  the  Hel 1 i nger-Rei ssner  principle  with  independent  strain. 
Initially  the  assumed  independent  strain  is  divided  into  a  lower  order  part  and 
a  higher  order  part.  The  stiffness  matrix  corresponding  to  the  lower  order 
assumed  strain  is  equivalent  to  the  stiffness  matrix  of  the  assumed  displacement 
model  element  with  the  reduced  integration  scheme.  The  spurious  kinematic  modes 
of  the  element  are  suppressed  by  introducing  a  stabilization  matrix  associated 
with  a  judiciously  chosen  set  of  higher  order  assumed  strain  fields.  Numerical 
results  show  that  this  element  is  free  of  locking  even  for  very  thin  plates  and 
shells. 


INTRODUCTION 

Since  the  early  days  in  the  history  of  the  finite  element  method,  a  great 
deal  of  research  effort  has  been  directed  to  the  finite  element  modeling  of  thin 
shell  structures.  Among  all  existing  approaches,  the  degenerate  solid  shell 
element  concept  [Ahmad,  Irons  and  Zienkiewicz  (1970)]  appears  to  be  the  most 
convenient  for  the  description  of  the  arbitrary  shell  geometry  and  the  kinema¬ 
tics  of  deformation.  However,  it  is  well  known  that  the  degenerate  solid  shell 
elements  exhibit  a  serious  drawback  unless  special  care  is  taken.  Tnis  phenome¬ 
non,  known  as  locking,  arises  from  the  overstiffening  effect  due  to  the  con¬ 
ditions  of  zero  inplane  strain  and  zero  transverse  shear  strain  when  the  shell 
thickness  becomes  small  [Lee  and  Pian  (1978)]. 

A  very  popular  way  of  alleviating  locking  has  been  to  utilize  the  reduced 
or  selective  integration  scheme  [Zienkiewicz,  Too  and  Taylor  (1971);  Pawsey  and 
Clough  (1971);  Hughes,  Cohen  and  Haroun  (1978);  Pugh,  Hinton  and  Zienkiewicz 
(1978);  Stolarski  and  Belytschko  (1982)].  However,  the  reduced  or  selective 
integration  scheme  has  not  been  successful  in  eliminating  the  effect  of  locking 


completely.  Even  with  the  2x2x2  point  reduced  integration,  an  eight  node  shell 
element  based  on  the  assumed  displacement  finite  element  model  still  experiences 
locking.  On  the  other  hand,  the  2x2x2  point  reduced  integration  rule  applied 
to  a  nine  node  element  or  the  3x3x2  point  reduced  integration  rule  applied  to 
a  sixteen  node  element  eliminates  the  effect  of  locking.  However,  they  intro¬ 
duce  spurious  kinematic  modes  which  lead  to  unstable  finite  element  models. 

To  improve  the  kinematic  stability,  we  may  employ  selective  integration  schemes 
to  these  elements  in  which  a  higher  order  integration  rule  is  used  for  the 
bending  part.  However,  selective  integration  schemes  cannot  eliminate  the 
unstable  spurious  kinematic  modes  completely.  In  short,  it  is  not  easy  to  find 
an  appropriate  reduced  or  selective  Integration  rule  which  can  eliminate  both 
locking  and  undesirable  kinematic  modes  at  the  same  time. 

In  order  to  suppress  the  spurious  kinematic  modes,  we  may  add  a  stabiliza¬ 
tion  matrix  to  the  element  stiffness  matrix  evaluated  by  a  reduced  integration 
rule  [Belytschko,  Ong  and  Liu  (1984);  Belytschko,  Liu,  Ong  and  Lam  (1985)].  In 
doing  so,  great  care  is  needed  to  avoid  reintroducing  the  effect  of  locking 
through  excessive  stabilization.  Recently,  a  rational  method  of  generating  a 
stabilization  matrix  has  been  developed  [Lee  and  Rhiu  (1986)].  This  method  is 
based  on  the  Hellinger-Reissner  principle  Including  both  independent  strain  and 
displacement-dependent  strain.  The  assumed  independent  strain  is  divided  into  a 
lower  order  part  and  a  higher  order  part.  With  a  proper  integration  rule,  the 
lower  order  assumed  strain  leads  to  an  element  stiffness  matrix  equivalent  to 
that  based  on  the  assumed  displacement  model  evaluated  with  the  same  integration 
rule  [Lee  (1978);  Malkus  and  Hughes  (1978)].  A  judiciously  chosen  higher  order 
independent  strain  field  is  used  to  generate  a  stabilization  matrix.  Following 
this  approach,  a  nine  node  element  which  is  free  of  locking  and  undesirabe 
spurious  kinematic  modes  has  been  developed  for  the  analysis  of  thin  shell 


structures  [ Rh i u  and  Lee  (1987);  Rhiu  (1985)]. 

Encouraged  by  this  success,  we  extend  in  this  paper  the  new  approach  to  the 
formulation  of  a  sixteen  node  degenerate  solid  shell  element.  Since  displacement 
fields  are  assumed  bicubic,  the  sixteen  node  element  has  the  potential  to  repre¬ 
sent  shell  behavior  with  considerable  accuracy.  However,  with  the  4x4x2  point 
integration,  the  sixteen  node  element  based  on  the  assumed  displacement  model 
still  suffers  from  locking,  particularly  for  distorted  or  curved  finite  element 
meshes.  On  the  other  hand,  as  mentioned  previously,  the  element  stiffness 
matrix  evaluated  by  the  3x3x2  point  reduced  integration  rule  has  unstable 
spurious  kinematic  modes.  These  spurious  kinematic  modes  will  be  identified. 

Then  they  will  be  suppressed  by  adding  a  stabilization  matrix  which  is  derived 
through  the  use  of  appropriately  assumed  higher  order  independent  strain  fields. 
Finally,  the  performance  of  the  present  element  will  be  tested  by  solving 
example  problems. 


GEOMETRY  AND  KINEMATICS 

Figure  1  shows  the  midsurface  of  a  curved  sixteen  node  shell  element.  In 
order  to  describe  the  shell  geometry  and  the  kinematics  of  deformation,  local 
coordinates  with  components  x,  y  and  z  are  defined  on  the  shell  midsurface  in 
addition  to  global  coordinates  with  components  X,  Y  and  Z.  The  x,  y  and  z  axes 
of  the  local  coordinate  system  are  parallel  to  the  orthogonal  unit  vectors  a^, 

$2  and  a^  respectively.  The  unit  vectors  a^  and  ^  are  tangential  to  the  shell 
midsurface  while  a,  is  normal  to  the  surface.  The  a,,  a0  and  a-j  vectors  are 
given  at  each  node  as  an  input.  In  addition,  they  are  defined  at  each  integra¬ 
tion  point  in  a  manner  which  will  be  discussed  later. 

With  the  coordinate  systems  described  above,  the  global  position  vector  X  of 
a  generic  material  point  can  be  expressed  as 

(1) 


3 


where  XQ  is  the  global  position  vector  of  a  point  located  on  the  shell  midsur- 
face,  5  j  £3  is  a  vector  drawn  from  the  point  on  the  midsurface  to  the  generic 
material  point,  t  is  the  shell  thickness  and  the  nondimensional  coordinate  c 
runs  from  -1  to  1.  Assuming  the  shell  undergoes  small  deformation,  the  displa¬ 
cement  vector  U  of  the  generic  material  point  with  respect  to  the  global  coor¬ 
dinate  system  can  be  expressed  as 

~*~o  +  57~~  (?) 


where 


&  =  l-ft2»a1l 


In  Eq.  (2b),  01  and  62  represent  small  rotations  of  a3  around  the  x  and  y  axes 
respectively.  In  Eq.  (2),  the  global  displacement  vector  UQ  of  the  point  on  the 
shell  midsurface  is  related  to  the  corresponding  local  displacement  vector  u 
with  components  u,  v  and  w  through  a  transformation  matrix  T  such  that 

Ho  *  X  SI  (3a) 

I  *  [il.i2.a3]  (3t>) 

Then  introducing  the  isoparametric  representation,  Eqs.  (1)  and  (2)  can  be 
expressed  as 


i  ,  IU  i 

X  -  l  N^c.n)  *1+^  N.(c,n)ti  a]  (4) 

16  .  16 

U  "  J  Mc.njT^  +  \  C  Ni(C.n)ti  ^  6.  (5) 

where  x’,  tit  53,  Xi ,  ,  bi ,  ei  are  the  values  of  XQ,  t,  43,  J,  jy .  b,  e  at  node 

1,  and  is  the  bicubic  shape  function  in  parent  coordinates  £  and  n. 

With  the  description  of  X  and  J.J  in  Eqs.  (4)  and  (5),  the  displacement- 


dependent  strain  vector  defined  with  respect  to  the  global  coordinate  system  can 
be  expressed  in  terms  of  the  vector  of  nodal  degrees  of  freedom.  Then,  using 
strain  transformation,  the  strain  vector  F  in  the  local  coordinate  system  is 
written  symbolically  as 

~  *  ^xx  ^yy  ^xy  ^yz  ^zx-^  ^ 

*  B(C,n,c) 

where  8(?,n,c)  is  the  strain-displacement  transformation  matrix  and  the 
element  nodal  degrees  of  freedom  vector  jge  Is  expressed  as: 

3e  3  Lj4i  ~1*  ~2  ~2 . .  ~16-J 

FINITE  ELEMENT  FORMULATION 

For  the  generation  of  our  sixteen  node  shell  element,  we  utilize  the 
Hel 1 inger-Reissner  functional  expressed  as  follows: 

*„  -  /  <et  c  r  -  4  et  c  E)  dV  -  U  (8) 

where  F  is  the  displacement-dependent  local  strain  vector  given  in  Eq.  (6)  and 
E  is  the  independent  local  strain  vector  such  that 

~  =  ^xx  ^yy  ^xy  ^yz  ^zx^  ^ 

In  the  present  formulation,  the  independent  strain  components  are  assumed 
to  be  linear  at  most  through  shell  thickness.  In  addition,  in  Eq.  (8),  W 
represents  the  applied  load  term,  V  is  the  volume  of  shell  and  C  is  a  5x5 
elastic  coefficient  matrix. 

Following  Lee  and  Rhlu  (1986),  initially  the  Independent  strain  E  is  divided 
into  two  parts  such  that 


5 


(10) 


~  s  ~L  + 

where  E^  is  the  independent  strain  vector  with  lower  order  assumed  polynomial 

terms  in  n  and  EH  is  the  higher  order  independent  strain  vector. 

Substituting  this  expression  into  Eq.  (8),  the  functional  *  becomes 

K 

*R  =  l  ue  -  w  (11) 

where 

ue  =  £  £  -  7  C  El)  dV  -  \  J  Ej  C  Eh  dV 

+  /  Ej  £  (I  -  £l)  dV  (12) 

and  l  indicates  summation  or  assembly  over  all  elements. 


For  a  sixteen  node  element  of  flat  rectangular  geometry,  the  displacement- 
dependent  strain  F  is  cubic  at  most  in  c  and  n.  If  the  lower  order  independent 
strain  £L  is  assumed  to  be  biquadratic  in  £  and  n,  the  first  integrals  in  Eq. 
(12)  can  be  integrated  exactly  in  £-n  plane  by  the  3x3  point  Gaussian  integra¬ 
tion  rule.  The  remaining  terms  are  integrated  by  the  4x4  point  rule  over  £  and 
n.  Although  these  integration  rules  are  determined  based  on  the  flat  rec¬ 
tangular  element  geometry,  the  same  integration  rules  will  be  adopted  for  ele¬ 
ments  with  arbitrary  geometry.  In  c-direction,  the  two  point  integration  rule 
is  used.  In  addition,  the  assumed  lower  order  independent  strain  can  be 
expressed  such  that 

18  18  _ 

■  j  NjU.n.c)  I1  -  I  N^C.n.t)  B(t i  ,n  1  ,c  1  )fle 


with 


=  BU.n.c)  ae 
18  _ 

B  -  I  N^c.n.c)  B(Ci,ni,Ci) 


(13) 


(13a) 


In  Eq.  (13),  shape  function  is  biquadratic  in  £,  n  and  linear  in  c  such  that 
N.  »  1  at  point  i  of  the  3x3x2  lower  order  integration  points  and  zero  at  other 


points,  and  E^  is  the  value  of  E  at  lower  order  integration  point  i.  Then,  for 
the  lower  order  strains,  it  is  possible  to  set 

iL  -  E  (14) 

at  the  3x3x2  integration  points. 

Applying  the  adopted  integration  rules  and  introducing  the  equivalence 
given  in  Eq.  (14),  Ug  in  Eq.  (12)  can  be  written  as 

ue  ■  1  {  il  £  k  dv  -  T  /H  si  £  Sh  m 

*  /  si  £  (I  -  SL)  dv  (15) 

H 

In  the  above  expression,  letters  L  and  H  under  the  integral  signs  represent  the 
lower  order  integration  (3x3x2  points)  and  the  higher  order  integration  (4x4x2 
points)  rules,  respectively. 

Based  on  the  limitation  principle  [Fraejis  de  Veubeke  (1965)],  the  polyno¬ 
mial  terms  in  the  assumed  strain  E^  cannot  be  of  higher  order  than  cubic  in  £ 
and  n.  Then,  with  biquadratic  JE^,  the  term  containing  in  the  last  integral 
of  Eq.  (15)  can  be  integrated  by  the  3x3x2  point  integration  rule.  Noting  this, 
Ug  can  be  rewritten  as 


Rhiu  and  Lee  (1987)  developed  a  nine  node  shell  element  using  the  expression 
for  Ug  equivalent  to  Eq.  (16).  For  the  present  sixteen  node  shell  element,  the 
expressions  In  Eq.  (15)  is  used. 

On  the  other  hand,  the  higher  order  assumed  strain  is  expressed  as 


7 


£h  “  £(«.**.«>  & 


where  £  is  the  assumed  strain  shape  function  matrix  which  contain  higher  order 
terms  in  £,  n  and  a  is  the  vector  of  higher  order  strain  parameters.  Note  that 


the  P  matrix  is  linear  in  c. 


Introducing  Eqs.  (6),  (13)  and  (18)  into  Eq.  (15),  the  functional  *R  in 


Eq.  (11)  becomes 


*R  =  I  (7  Se  ~L  3e  +  ~  £  3e  ”  \  ~  ~  ~  "  3e  iM 


where 


K.  =  /  BTC  B  dV 
~L  '  ^  ~  ^  ~ 


G  «  J  PT  C  (B  -  ¥)  dV 

^  #y  /V  IV  IV 

H 


H  -  /  PT  C  P  dV 
H  ~  ~  ~ 


l  fll  QP  =  w 


Setting  6 =  0  with  respect  to  a  results  in  the  compatibi 1 ity  equation  in 


discretized  form  as  follows: 


a  «  H  G  q. 


for  each  element. 


By  introducing  Eq.  (24)  into  Eq.  (19),  *R  can  be  written  as 


=  l  (7  3e  ~e  5e  "  3e  9e) 


In  the  above  equation,  the  element  stiffness  matrix  K  is  given  as 


*e  =  *L  +  *s 


where 


<  •  sT  tf 1  fi 


The  matrix  is  evaluated  by  the  3x3x2  point  integration  rule  while  the 
matrix  associated  with  the  higher  order  assumed  strain  is  evaluated  by  the  4x4x2 


point  integration  rule.  Note  that  the  matrix  is  in  fact  the  same  element 
stiffness  matrix  derived  from  the  conventional  assumed  displacement  model  based 
on  the  principle  of  virtual  work  with  the  3x3x2  point  reduced  integration  rule. 
The  Kl  matrix  has  spurious  kinematic  modes,  and  these  modes  are  suppressed 
by  adding  a  properly  constructed  <s  matrix.  Thus,  K$  plays  the  role  of  a  sta¬ 
bilization  matrix. 

To  construct  the  element  stiffness  matrix,  it  is  necessary  to  evaluate  the  B 
matrix  at  both  the  higher  order  integration  points  and  the  lower  order  integra¬ 
tion  points.  Alternately  B  at  the  lower  order  integration  points  can  be  inter¬ 
polated  from  B  evaluated  at  the  higher  order  integration  points  as  follows: 


~  ^i »  ni  »ci ) 


32 

=  I 

j  =  l 


*n  ^  •C'j) 


8Uj,nj,Cj) 


(28) 


where  the  subscripts  i  and  j  stand  for  the  lower  order  integration  points  and 
the  higher  order  integration  points  respectively,  and  the  shape  function  N.  is 

J 

bicubic  in  c,  n  and  linear  in  c  such  that  N.  =  l  at  the  point  j  of  the  4x4x2 

J 

higher  order  integration  points  and  zero  at  other  points. 


CONTROL  OF  THE  SPURIOUS  KINEMATIC  HOPES 
For  an  element  of  flat  rectangular  shape  with  sides  along  x  *  ±1  and  y  =  ±1 
lines.  It  Is  possible  to  determine  the  analytical  expressions  for  the  spurious 
kinematic  modes  of  the  matrix  by  expressing  the  assumed  u,  v,  w,  0^  and  02  as 
polynomial  functions  in  x  and  y  coordinates.  For  example,  we  may  write 

u  =  a:  +  a2x  +  .  +  a16x3y3  (29) 

and  similarly  for  v,  w,  0^  and  02«  Then  displacement-dependent  strain  vector 
E  can  be  expressed  from  these  assumed  displacement  fields.  Now  noting  that 
spurious  kinematic  modes  of  the  K,  matrix  do  not  produce  strain,  we  set 


(30) 


E  =  0 

at  the  3x3x2  lower  order  integration  points.  This  leads  to  a  set  of  72  homoge¬ 
neous  equations  from  which  we  can  identify  the  following  seven  spurious  kinema¬ 
tic  modes: 


(1)  u  =  -Cj  y  (3x2  -  5x2y2  +  y2 ) 
v  =  x  (x2  -  5x2y2  +  3y2) 

(2)  ei  -  c5  x  (-  f  +  x2  -  5x2y2  +  3y2) 

02  =  c5  y  (■  t  +  3x2  *  5*V  +  y2) 

(3)  u  =  C2  xy  (9  -  15x2  -  15y2  +  25x2y2) 

(4)  v  =  C3  xy  (9  -  15x2  -  15y2  +  25x2y2 ) 

(5)  w  =  C4  xy  (9  -  15x2  -  15y2  +  25x2y2 ) 

(6)  6j  *  C6  xy  (9  -  15x2  -  15y2  +  25x2y2 ) 

(7)  *  c7  *y  (9  -  15x2  *  1 5y2  +  25x2y2) 


(31a) 

(310) 

(31c) 
( 3  Id ) 
( 31e) 
( 31  f ) 
(31g) 


where  Cj,  C2,  .....  C7  are  arbitrary  constants.  The  modes  given  in  Eqs.  (31a) 
and  (31b)  are  incompatible.  That  is,  they  disappear  for  an  assembly  of  only  two 
elements.  However,  the  remaining  modes  are  compatible  and  persist  even  after 
assembling  elements,  resulting  in  an  unstable  finite  element  model.  These 
spurious  kinematic  modes  are  suppressed  by  introducing  carefully  chosen  higher 
order  assumed  strain  fields  as  follows: 

The  displacement-dependent  strain  component  corresponding  to  Eqs.  (31a)  to 
(31g)  are 


Exx  =  -  Cj  (6xy  -  10xy3)  +  C5  z(6xy  -  10xy3 ) 
+  C2  (9y  -  45x2y  -  15y3  +  75x2y3) 

+  C?  z(9y  -  45x2y  -  15y3  +  75x2y3) 


(32a) 


Eyy  =  C1  ^6xy  "  10x  y*  ‘  C5  z^6xy  "10x  y* 

+  C3  (9x  -  15x3  -  45xy2  +  75x3y2) 

-  Cg  z(9x  -  15x3  -  45xy2  +  75x3y2 )  (32b) 

Eyv  -  C9  (9x  -  15x3  -  45xy2  +  75x3y2) 

Ajr  c 

+  C3  (9y  -  45x2y  -  15y3  +  75x2y3 ) 

+  C7  z(9x  -  15x3  -  45xy2  +  75x3y2) 

-  Cg  z(9y  -  45x2y  -  15y3  +  75x2y3)  (32c) 

EyZ  =  C4  (9x  -  15x2  -  45xy2  +  75x3y2) 

-  C5  (-  |  x  +  x3  -  5x3y2  +  3xy2 ) 

-  Cg  (9xy  -  15x3y  -  15xy3  +  25x3y3)  (32d) 

Ezx  *  C4  (9y  -  45x2y  -  15y3  +  75x2y3 ) 

♦  C5  (-f  y  ♦  3x2y  -  5x2y3  ♦  y3 ) 

+  C?  (9xy  -  15x3y  -  15xy3  +  25x3y3)  (32e) 

Examining  Eqs.  (32a)  to  (32e),  we  realize  that  the  spurious  kinematic  modes  In 
Eqs.  Ola)  to  (31g)  are  suppressed  for  the  following  higher  order  assumed  strain 
fields: 

ttxx)H  a  °i  x  y  +  a5  y  +  «7  *y  +  <*8  zxy 
(Eyy>H  ■  “2  *  °6  lxV*  °9  *  “10  zx’» 


!e,Jh  *  o 


(33) 


In  Eq.  (33),  o^,  «2,  ••••  ••  a1Q  are  unknown  coefficients.  Alternately,  noting 
that  the  modes  corresponding  to  and  Cg  are  incompatible,  we  may  drop 
Og,  ag  and  a^g  terms  from  Eq.  (33).  This  leads  to  an  assumed  higher  order 
strain  field  with  six  coefficients  and  the  resulting  element  stiffness  matrix 
has  eight  zero  eigenvalues.  However,  when  elements  are  assembled,  the  resulting 
finite  element  model  is  kinematically  stable. 

For  an  element  with  arbitrary  geometry,  we  use  5  n  ,  5  n  terms  etc. 

2332  23  32 

instead  of  x  y  ,  x  y  terms  etc.  Since  £  n  and  £  n  are  not  symmetric  with 

respect  to  parent  coordinates  £  and  n,  the  element  stiffness  matrix  may  be 
dependent  on  the  choice  of  local  coordinate  systems  used.  If  the  local  coor¬ 
dinate  system  is  chosen  such  that  the  a1  or  x  axis  is  parallel  with  the  £  coor¬ 
dinate,  then  the  element  stiffness  matrix  is  not  invariant  when  element  geometry 
is  nonrectangular.  For  example,  consider  the  distorted  elements  with  different 
node  numberings  as  shown  in  Fig.  2.  Even  though  both  elements  have  the  same 
geometry,  we  obtain  two  different  element  stiffness  matrices.  The  local  coor¬ 
dinate  system  with  x  or  parallel  with  £  has  been  used  in  Lee,  Wong  and  Rhiu 
(1985)  in  conjunction  with  a  nine  node  shell  element.  In  spite  of  the  lack  of 
invariance  of  the  element  stiffness  matrices,  this  nine  node  shell  element 
showed  excellent  performance.  This  indicates  that  the  invariance  property  is 
not  absolutely  necessary  for  a  good  finite  element.  However,  in  the  present 
study,  we  enforce  the  invariance  of  element  stiffness  matrices  by  assigning  a 
particular  local  coordinate  system  for  a  given  geometry  of  elment  as  follows 
[Rhiu  and  Lee  (1987)] : 

If  XQ  denotes  the  position  vector  of  the  point  located  at  £  *  n  *  c  =  0,  we 
may  define  two  unit  vectros  v^  and  y2  at  this  point  such  that 


The  angle  0Q  between  these  two  unit  vectors  is  determined  by  the  following 
equation, 

cos  8  =  v,  •  v0 

o  ~i  ~z 


(35) 


Then,  if  eQ  is  less  than  or  equal  to  90°,  the  unit  vector  a^  in  the  x  direction 
of  local  coordinate  system  is  chosen  to  be  parallel  to  £  axis  such  that 


(36a) 


Otherwise,  Sj  is  parellel  to  n  axis  such  that 


(36b) 


With  this  choice  of  a^,  we  can  easily  determine  the  other  two  unit  vectors 
&2  and  43,  with  ^  being  normal  to  the  shell  midsurface.  Note  that,  while 

and  V£  are  determined  at£*n*cs0  point,  the  a^,  £2  and  a-j  vectors  can  be 
computed  at  any  point  on  the  shell  midsurface.  In  particular,  a^,  ^  and  a-j  are 
needed  at  the  integration  points. 

With  the  local  coordinate  system  defined  as  above,  the  higher  order  assumed 
strains  for  the  sixteen  node  shell  element  are  chosen  as  follows: 


£h 


*  P  o 


(37) 


where  for  the  10a  version. 


0  9j  0  ;gj  0  0 


g2  0 


P  =  0  0  0  0  0  0 


0  0  0  0  g1  0 


0  0  0  0  0  f. 


0  0 


i 

“  =  Lai»  <*2 . ai&i 


and  for  the  6a  version. 


^  0  cfj  0  0  0 


o  9, 


0  0 


P  «  0  0 


0  0  0 


0  0 


o  9j  0 


0  0 


0  0  f. 


i 

L°l»  a2*  . “(jJ 


In  Eqs.  (38a)  and  (39a),  f^,  g^,  and  g^  are  chosen  as  follows: 

(1)  if  x  or  a^  is  parallel  to  e  as  in  Eq.  (36a) 


,  ,.2  3-  .  3 

f,  a  C  n  ,  fo  ■  £n 


3  2  3 

9j  =  e  n  ,  g2  -  C  n 


(2)  if  x  or  aj  is  parallel  to  n  as  in  Eq.  (36b) 


,  .3  2  -  „3 

fi  *  M  ,  f,  *  5  1 


9j  *  C2n3,  g2  -  Cn3 


NUMERICAL  TESTS 


In  order  to  evaluate  the  performance  of  the  present  sixteen  node  element. 


several  numerical  tests  involving  simple  plates  and  shells  were  carried  out. 


For  the  purpose  of  identification,  the  present  sixteen  node  element  is  called 
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SHEL16.  Whenever  possible,  the  effectiveness  of  SHEL16  element  is  compared  with 
the  DISP16  element  based  on  the  conventional  displacement  model  with  the  4x4x2 
point  integration  rule.  Host  of  the  numerical  results  are  presented  in  tabular 
form  so  that  they  can  be  used  for  future  reference.  For  the  SHEL16  element,  it 
turns  out  that  numerical  results  for  the  6a  assumed  strain  and  the  10a  assumed 
strain  are  almost  identical  for  the  cases  tested  in  this  paper.  Therefore  only 
the  results  for  the  6a  version  is  presented.  All  numerical  examples  were 
calculated  with  double  precision  accuracy  on  the  UNIVAC  1100/92  machine  at  the 
University  of  Maryland. . 

(a)  A  Simply  Supported  or  Clamped  Square  Plate 

Plate  bending  problems  provide  examples  to  investigate  the  effect  of  trans¬ 
verse  shear  locking  alone.  A  quarter  of  a  square  plate  subjected  to  uniformly 
distributed  load  p  was  modeled  by  uniform  lxl  and  2x2  meshes  and  distorted  2x2 
and  4x4  meshes  as  shown  in  Figs.  3(a)  to  3(c).  Both  simply  supported  and 
clamped  boundary  conditions  were  considered. 

Table  1  lists  the  computed  nondimensional  deflection  at  the  centroid  of  the 
plate.  These  values  are  normalized  with  respect  to  the  analytical  solution 
based  on  the  Kirchhoff  thin  plate  theory  [Timoshenko  and  Woi nowski-Krieger 
(1959)].  For  the  simply  supported  plate,  both  SHEL16  and  DISPlo  elements  give 
numerical  results  very  close  to  the  analytical  solutions  for  the  uniform  meshes. 
For  the  distorted  2x2  mesh,  the  SHEL16  element  does  not  suffer  any  transverse 
shear  locking  over  a  wide  range  of  L/t  ratios  while  the  0ISP16  element  reveals  a 
slight  effect  of  shear  locking  when  the  plate  becomes  very  thin.  For  the 
distorted  4x4  mesh,  both  elements  give  very  accurate  results.  For  the  clamped 
plate,  the  SHEL16  element  gives  very  accurate  and  reliable  numerical  results 
over  a  wide  range  of  L/t  ratios  regardless  of  mesh  distortion.  However, 
for  the  distorted  meshes,  the  performance  of  the  DISP16  element  deteriorates  as 


the  plate  becomes  thin.  Even  in  this  case,  the  4x4  mesh  shows  very  accurate 
solution  up  to  L/t  =  10,000. 

2 

Table  2  shows  nondimensional  bending  moments  Mx/pL  per  unit  length  eva¬ 
luated  at  integration  point  E  and  nondimensional  shear  forces  Qx/pL  per  unit 
length  evaluated  at  integration  point  F.  Note  that  the  SHEL16  element  solu¬ 
tions  are  totally  insensitive  to  the  wide  range  of  L/t  ratios  considered  here. 
Table  2  also  includes  analytical  solutions  obtained  at  corner  points  C  and  D. 
They  are  listed  to  check  the  order  of  magnitude  of  numerical  solutions. 

(b)  A  Pinched  Cylindrical  Shell 

As  a  deep  shell  example,  a  cylindrical  shell  loaded  at  two  opposite  poincs 
as  shown  in  Fig.  4  was  tested.  Both  diaphragmed  and  fixed  edge  conditions  were 
considered.  Due  to  symmetry  in  geometry  and  loading,  only  one  octant  of  the 
shell  was  modeled  by  3x4,  4x5  and  5x6  meshes  as  shown  in  Figs.  5(a)-(c).  In 
addition,  as  shown  in  Fig.  5(d),  an  irregular  mesh  designated  as  5x61  was  also 
considered.  Note  that  the  meshes  illustrated  in  Figs.  5(a)-(d)  are  on  the 
stretched  plane  of  the  octant  ABCO  of  the  shell.  Moreover,  in  order  to  describe 
more  accurately  the  complex  shell  behavior  In  the  region  near  the  load  point  C, 
fine  meshes  are  used  along  lines  BC  and  CD. 

Table  3  lists  the  nondimensional  displacements  at  various  points  on  the 
diaphragmed  shell  for  R/t  =  100,  300  and  500.  They  are  compared  with  the  analy¬ 
tical  solutions  given  by  Fliigge  (1962).  The  analytical  solution  is  based  on  a 
shell  theory  which  neglects  the  effect  of  transverse  shear  deformation.  Table  3 
also  includes  numerical  results  obtained  by  the  DISP16  element  with  the  5x6 
mesh.  For  the  models  with  SHEL16  elements,  the  solutions  get  closer  to  the  ana¬ 
lytical  solutions  as  the  number  of  elements  increases.  It  is  noteworthy  that 
the  solutions  for  the  distorted  5x61  mesh  are  very  close  to  that  for  the  regular 
5x6  mesh.  On  the  other  hand,  the  DISP16  element  shows  signs  of  locking  as  the 
solutions  deteriorate  with  increasing  R/t  ratios.  Even  for  R/t  *  100,  the  5x6 


mesh  solution  with  the  DISP16  element  is  worse  than  the  3x4  mesh  solution  witn 
the  SHEL16  element. 

Table  4  lists  nondimensional  deflections  at  the  pinched  point  C  of  the  shell 
with  fixed  ends.  A  good  convergence  is  observed  as  the  finite  element  model 
with  the  SHEL16  elements  is  refined.  Also  there  is  no  significant  discrepancy 
between  the  5x6  mesh  and  the  5x61  mesh.  Figs.  6  and  7  show  inplane  force 
and  moment  per  unit  length  along  line  BC  for  the  5x6  mesh  with  SHEL16  ele¬ 
ments.  An  analytical  solution  for  the  fixed  ends  case  is  not  available. 

(c)  A  Hemispherical  Shel  1 

As  a  doubly  curved  shell  example,  a  hemispherical  shell  subjected  to  con¬ 
centrated  loads  as  shown  in  Fig.  8  was  considered.  This  problem  exhibits  pre¬ 
dominantly  bending  behavior  with  very  little  inplane  behavior.  Due  to  symmetry 
in  geometry  and  loading,  a  quarter  of  the  shell  was  modeled  by  4  element,  9  ele¬ 
ment,  16  element  and  20  element  meshes.  The  4  element,  9  element  and  16  element 
meshes  are  created  by  dividing  uniformly  over  the  angles  6  and  t.  The  20  ele¬ 
ment  mesh  is  created  from  the  16  element  mesh  as  shown  in  Fig.  8. 

For  convenience,  a  small  region  at  point  C  was  not  Included  in  the  finite 
element  modeling.  As  a  check,  two  different  cases  were  tested.  In  one  case,  the 
region  within  0  *  0.5°  was  cut  out  while,  in  the  other  case,  the  region  within  e 
=  1°  was  excluded.  The  two  cases  gave  the  same  result.  In  table  5  the  com¬ 
puted  nondimensional  deflection  DW^/PR2  at  point  A  is  compared  with  the  analyti¬ 
cal  solution  reported  by  Morley  and  Morris  (1978).  Symbol  D  represents  bending 
rigidity.  The  analytical  solution  is  based  on  the  Rayleigh-Ritz  method.  For 
R/t  *  250,  the  solution  for  the  16  element  model  agrees  exactly  with  the  analy¬ 
tical  value  of  0.185.  Even  the  10  element  model  shows  only  0.05%  error.  On  the 
other  hand,  the  DISP16  element  suffers  from  locking.  Table  5  also  includes  the 
R/t  *  500  case.  Morley  and  Morris  (1978)  did  not  consider  this  case. 


(d)  A  Toroidal  Shell  under  Internal  Pressure 

A  toroidal  shell  subjected  to  an  internal  pressure  p  was  analyzed  by  the 
SHEL16  element.  Figures  9(a)  and  9(b)  show  the  geometry  and  material  data.  Tne 
toroidal  shell  has  both  positive  and  negative  curvatures  along  the  meridional 
angle.  Due  to  the  horizontal  plane  of  symmetry  and  the  axisymmetric  loading,  an 
upper  sector  of  shell  with  an  angle  of  8°  was  modeled  with  a  row  of  13  elements 
and  a  row  of  26  elements.  The  subtended  angles  of  individual  elements  in  the 
13  element  model  are  listed  in  Table  6.  The  26  element  model  is  obtained  from 
the  13  element  model  by  dividing  each  element  into  two  elements  with  equal  sub¬ 
tended  angles.  Numerical  results  for  the  13  element  model  and  the  26  element 
model  were  almost  the  same.  Therefore  only  the  13  element  solutions  are 
reported  here. 

Table  7  shows  nondimensional  normal  deflection  (w/r)  x  10  in  comparison 
with  the  numerical  solution  by  Kalnins  (1964)  for  r/t  ratios  of  20  and  200. 
Kalnins'  solution  is  a  combination  of  the  direct  integration  and  the  finite  dif¬ 
ference  method.  A  very  good  agreement  between  the  results  of  the  SHEL16  element 
and  Kalnins'  solution  is  observed.  Table  7  also  includes  the  SHEL16  element 
solution  for  r/t  *  1,000.  This  case  was  not  considered  by  Kalnins.  Figures  10 
and  11  show  the  deflection  and  the  bending  stress  (o0e)b  at  the  top  surface  of 
the  shell  along  the  meridional  angle  direction.  An  excellent  agreement  between 
the  two  solutions  is  observed  for  r/t  =  20  and  200.  For  r/t  =  200,  the  distri¬ 
bution  of  the  membrane  stress  (oe0)m  is  shown  in  Fig.  12.  Again  the  SHEL16  ele¬ 
ment  solution  is  almost  identical  to  Kalnins'  solution.  The  (o„.)  /E  curves  for 
r/t  =  20  and  1,000  are  very  close  to  that  for  r/t  *  200.  Therefore,  they  are 
not  shown  in  Fig.  12  to  avoid  cluttering. 

CONCLUSION 

Results  of  numerical  tests  demonstrate  that  the  present  SHEL16  element 
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can  be  used  to  provide  reliable  solutions  for  thin  plates  and  shells  regardless 
of  distorted  element  geometries  and  clamped  boundary  conditions.  In  addition, 
the  SHEL16  element  with  the  10a  version  assumed  strain  is  kinematically  stable 
at  element  level  while  the  SHEL16  element  with  the  6a  version  assumed  strain  is 
kinematically  unstable  at  element  level  but  stable  at  global  structural  level. 
Thus  the  stabilization  scheme  with  a  judiciously  chosen  set  of  higher  order 
assumed  strain  terms  has  successfully  suppressed  compatible  kinematic  modes 
without  reintroducing  the  locking  effect.  Finally,  the  present  SHEL16  element 
can  be  used  to  generate  benchmark  solutions  for  testing  the  performance  of  other 
shell  elements. 
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Table  1, 


Maximum  nondimensional  deflection  at  the  centroid  of  the 
square  plate  under  uniform  pressure 


PI  ate 


Uni  form: 
lxl 


Simply 


Supported 


Distorted: 

2x2 


Uni  form: 
lxl 


Clamped 


Di stored: 
2x2 


Element 


SHEL16 

1.0000 

0.9995 

0.9995 

0.9995 

DISP16 

1.0153 

1.0150 

1.0150 

1.0150 

SHEL16 

1.0005 

1.0000 

1.0000 

1.0000 

DISP16 

1.0012 

1.0007 

1.0007 

1.0007 

SHEL16 

1.0007 

1.0007 

1.0005 

1.0005 

0ISP16 

1.0015 

0.9956 

0.9542 

0.9380 

SHEL16 

1.0005 

1.0002 

1.0002 

1.0002 

DISP16 

1.0000 

1.0002 

1.0000 

0.9906 

SHEL16 

0.9968 

0.9945 

0.9945 

0.9945 

DISP16 

1.0482 

1.0474 

1.0474 

1.0474 

SHEL16 

1.0024 

1.0000 

1.0000 

1.0000 

DISP16 

1.0016 

1.0000 

1.0000 

1.0000 

SHEL16 

1.0024 

1.0000 

0.9960 

0.9960 

DISP16 

0.9945 

0.9486 

0.3007 

0.0048 

SHEL16 

1.0024 

1.0000 

1.0000 

0.9992 

DISP16 

1.0016 

0.9992 

0.9929 

0.8632 

Table  2 


Nondimensional  bending  moment  and  shear  force  of  the 
square  plate  {SHEL16  element  with  uniform  2x2  mesh) 


Simply  Supported  Plate 

Clamped 

Plate 

M 

Q 

M 

Q 

L 

t 

(—) 

PL2  E 

(pU}F 

V,E 

(pf}l 

102 

0.0476 

0.310 

0.0226 

-0.41 

103 

0.0476 

0.310 

0.0226 

-0.41 

io4 

0.0476 

0.310 

0.0226 

1 

o 

• 

io5 

0.0476 

0.310 

0.0226 

• 

o 

i 

M  Q  M 

Analytical  Mj>c  (p£>0  (^C 


0.0479 


0.338 


0.0231 


Table  3.  Nondimensional  displacements  for  the 
pinched  cylinder  with  diaphragmed  ends 


R 

t 

Element 

Mesh 

_  Et"c 
p 

EtwB 

P 

.  Et“D 

P 

3x4 

165.3 

0.6776 

4.102 

SHEL16 

4x5 

166.1 

0.5218 

4.113 

5x6 

166.3 

0.4770 

4.113 

5x61 

166.3 

0.4718 

4.113 

100 

DISP16 

5x6 

159.1 

1.497 

4.087 

Analytical 

164.3 

0.4693 

4.114 

3x4 

636.3 

12.52 

9.778 

SHEL16 

4x5 

642.3 

12.52 

9.785 

5x6 

646.9 

12.22 

9.853 

5x61 

646.5 

12.32 

9.853 

300 

0ISP16 

5x6 

531.1 

21.31 

9.397 

Analytical 

647.3 

9.867 

3x4 

1172.2 

10.60 

14.46 

SHEL16 

4x5 

1200.9 

10.23 

14.45 

5x6 

1212.0 

13.47 

14.55 

5x61 

1210.2 

13.21 

14.54 

500 

DISP16 

5x6 

847.6 

5.904 

13.32 

Analytical 

1223.4 

14.67 
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Table  6.  Meridional  subtended  angle  of  the  elements  for  the  toroidal  shell 


element 


no. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

A0 

20 

20 

20 

11 

10 

7 

4 

7 

10 

11 

20 

20 

20 

(degrees) 


m 


S 


Nondimensional  deflections  (w/r)  x  10  of  the  toroidal  shell 


r/t 

=  20 

r/t 

«  200 

r/t  *  1000 

SHEL16 

Kalnins 

SHEL16 

Kalnins 

SHEL16 

0.1034 

0.103 

0.1038 

0.100 

0.1034 

4.2119 

4.208 

5.1423 

5.151 

5.2911 

3.4668 

3.467 

3.2952 

3.297 

2.8179 

1.2513 


1.249 


1.3038 


1.298 


1.3273 


E:  X  =  0.056in 
Y  =0.056in 


F:  X  =0.944in 
Y  =0.056  in 


E 


L/2 


L=  2in 

E  =  I07  psi  v  =  0.3 


Figure  3(a)  A  square  plate:  2x2  uniform  mesh 


-  R/t  =  100 

- R/t  =300 

- R/t  =  500 
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lure  6  Inplane  force  along  line  BC  of  the  pinched  cylindrical  shell 
with  fixed  ends 


Figure  7  Bending  moment  distribution  along  line  BC  of  the  pinched 
cylindrical  shell  with  fixed  ends 


A  X 


R/t  =250,500 
E  =  1.0,  v  =0.3 
P  =  2.0 


Figure  8  Finite  element  mesh  for  one  quarter  of  the  hemispherical  shell 
subjected  to  concentrated  loads 


toroidal  shell  under  internal  pressure 


SECTION 


toroidal  shell  under  internal  pressure  p:  section 


KALNINS  r/t  =20,200 
SHELI6  r/t  =20 


deflection  along  meridional  angle  direction  for  the  toroidal  shell 


KALNINS  r/t  =20,200 ;- 
SHELI6  r/t  =20 

=200  = 


g  stress  along  meridional  angle  direction  for  the  toroidal  shell 


